Resonance instability of axially-symmetric magnetostatic equilibria 
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We review the evidence for and against the possibility that a strong enough poloidal field stabilizes 
an axisymmetric magnetostatic field configuration. We show that there does exist a class of resonant 
MHD waves which produce instability for any value of the ratio of poloidal and toroidal field strength. 
We argue that recent investigations of the stability of mixed poloidal and toroidal field configurations 
based on 3-d numerical simulations, can miss this instability because of the very large azimuthal 
wave numbers involved and its resonant character. 

PACS numbers: 47.20.-k, 47.65.-d, 95.30.Qd 



I. INTRODUCTION 

The stability of hydromagnetic configurations is still 
a topic of debate. Even simple magnetic configura- 
tions consisting in a pure azimuthal (toroidal) or ver- 
tical (poloidal) field are generally unstable (see, e.g., 
yet the magnetic fields observed in several astrophysical 
contexts are stable on a secular time scale. In this 
context, the energy principle of Bernstein et al. @ has 
extensively been used in the past to study the stability of 
simple poloidal or toroidal fields 043 and also of mixed 
combinations of the two [6|. In cylindrical geometry, it 
can be proved that the plasma is stable for all azimuthal 
and vertical wave numbers to and k, if it is stable for 
to = in the k — > limit, and for to = 1 for all k [Tj. On 
the other hand, to show that a generic configuration with 
a combination of vertical field and non-homogenous az- 
imuthal field is stable against the to = 1 mode (for all k) 
is not an easy task in general and one has to resort either 
to a variational approach or to a numerical investigation 
of the full eigenvalue problem in the complex plane @. 
In this respect, the "normal mode" approach can be more 
useful in astrophysics, as it is often important to know 
the growth rate of the instability and the properties of 
the spectrum of the unstable modes [t| [l(| • 

In recent years, the use of 3D numerical simulations 
has opened up the possibility of studying the stability of 
various field configurations following the evolution from 
the linear phase to the non-linear regime. A strategy of- 
ten used is to evolve a generic initial state which eventu- 
ally relaxes to a final configuration assumed to be stable 
[ll| - [l5| . The drawback with this approach is that it is 
difficult to characterize the topology of the final config- 
uration from the analysis of the numerical data and to 
determine a class of sufficient conditions for instability 
which could be of astrophysical interest. In particular, 
the conclusions of some recent works in this direction 
seem to point out that it is the strength of the poloidal 
field which stabilizes the basic state [lj, [l4| . 

The aim of this paper is to clarify that field configu- 
rations containing generic combinations of axial and az- 
imuthal fields are subject to a class of resonant MHD 



waves which can never be stabilized for any value of the 
ratio of poloidal and toroidal fields. The instability of 
these waves has a mixed character, being both current- 
and pressure-driven (la ]. We argue that in this case the 
most dangerous unstable modes are resonant, i.e. the 
wave vector k — (m/s)eg + k z e z is perpendicular to the 
magnetic field, B-k = where k z is the wavevector in the 
axial direction, to is the azimuthal wavenumber, and s is 
the cylindrical radius. The length scale of this instabil- 
ity depends on the ratio of poloidal and azimuthal field 
components and it can be very short, while the width of 
the resonance turns out to be extremely narrow. For this 
reason its excitation in simulations can be problematic. 

The paper is organized as follows. In Sec. 2, the main 
equations governing the behaviour of linear perturbations 
in cylindrical plasma configurations are presented. In 
Sec. 3, we consider a linear stability analysis of such con- 
figurations, using an analytical approach complemented 
by a numerical investigations. Direct numerical simula- 
tions of the non-linear evolution of a cylindrical configu- 
ration are presented in Sec. 4. In Sec. 5, we compare our 
results with those obtained by other authors and discuss 
possible astrophysical applications of this instability. 



II. BASIC EQUATIONS 

Let us consider an axially symmetric basic state with 
azimuthal and axial magnetic fields. The azimuthal field 
is assumed to be dependent on s alone, B v = B v (s), but 
the axial magnetic field B z is constant. We assume that 
the sound speed is significantly greater than the Alfven 
velocity in order to justify the use of incompressible MHD 
equations 

dv .„ w VP (VxB)xB 
at p Airp 



dB 
~dt 



Vx(vxB) = 0, V-B = 0, V-tf=G. (1) 



In the basic state, hydrostatic equilibrium in the ra- 
dial direction is assumed. We study a linear stabil- 
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ity with respect to small disturbances. Since the basic 
state is stationary and axisymmctric, the dependence of 
disturbances on t, tp, and z can be taken in the form 
exp (at — ik z z — imip). Linearizing Eq.(l) and eliminat- 
ing all variables in favor of the radial velocity distur- 
bance, vis, we obtain 
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where ua — (B ■ k)/^/Wp, uja z = k z B z /y/A-Kp, lu b = 
B v / s^JAitp, a = 91n_B/91ns, and A = 1 + m 2 /s 2 k 2 . 
Eq.([2]) describes the stability problem as a nonlinear 
eigenvalue problem. This equation has been first derived 
by Freidberg [l?} in his study of MHD stability of a dif- 
fuse screw pinch (see also [l(| )■ The author found that, 
for a given value of k z , it is possible to obtain multiple 
values of the eigenvalue er, each one corresponding to a 
different eigenfunction, and calculated a for the fastest 
growing fundamental mode. The most general form of 
Eq.(2), taking into account compressibility of plasma, 
was derived by Goedbloed [18j. Since we study the sta- 
bility assuming that the magnetic energy is smaller than 
the thermal one, the incompressible form of Eq.Q can 
be a sufficiently accurate approximation. In fact, Eq.Q 
was studied by Bonanno & Urpin in their analysis of 
the non-axisymmetric stability of stellar magnetic fields. 

We can represent the azimuthal magnetic field as B v = 
B V Qip(s), where B v q is its characteristic strength and 
ip(s) ~ 1. It is convenient to introduce dimensionless co- 
ordinate x = s/s2 and dimensionless quantities q — k z S2, 
T = <t/ljbo, w B0 = B vQ /s 2 y / 4:Trp, and e = B z /B v0 . 
Then, Eq. (2) transforms into 
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where 



/ = (?£• 



ip(x) 



A = 



q 2 x 2 (T 2 +f 2 ) 



q 2 x 2 



(4) 



With appropriate boundary conditions, Eq. (3) allows 
to determine the eigenvalue T. If the inner boundary is 
extended to include the cylinder axis it is not difficult 
to show that the eigenfunction for m = 1 must be non- 
vanishing there to ensure regularity. This result follows 
from the series solution of Eq.(3) near x = 0, so that 
vis oc x b with b — — l±m, and regularity at x = implies 



b = for m = 1, and b > for m > 1. In the setup 
discussed in this paper the inner boundary is not located 
at the axis, and we can safely assume that Vi s = at 
x = x% and x — x%. We will demonstrate the occurrence 
of a resonance instability in magnetic configurations by 
an analytical and numerical solution of Eq. (3), and by 
3D direct numerical simulations. 



III. LINEAR ANALYSIS OF INSTABILITY 

A. Analytical considerations 

It is interesting to have a qualitative understanding of 
the MHD spectrum, thus solving Eq. (3) in the small 
gap approximation. In this case one assumes that the 
distance between the boundaries, Aa; = x% — x\, is small 
compared to Xi = 1 and neglect in Eq. (3) terms of the 
order vi s /x compared to dvi s /dx. In this approximation, 
all coefficients of Eq. (3) can be considered as constant 
and Eq. (3) yields 
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The solution, satisfying the boundary conditions, is vi s oc 
sin[7r(a; — xi)/ Ax]. The corresponding dispersion relation 
is biquadratic and can be easily solved. The solution is 
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where p = q 2 /[q 2 + m 2 + (71-/A2;) 2 ]. The parameter / 
characterizes departures from the magnetic resonance, 
lja = 0. To show the occurrence of instability, we con- 
sider solution (6) at small departures from the magnetic 
resonance, / w 0. If a > 1, we have 



2m 2 (a - 1) 
! + (p 2 + m 2 )e 2 



(7) 



where p 2 = (ir/Ax) 2 and the instability is never sup- 
pressed for any finite value of e. The growth rate is a 
rapidly increasing function of m and T 2 s» (1 + e 2 ) -1 in 
the limit m^> p 2 . If a < 1, then Eq. (6) yields 



r 



a 



1 — a ' 



(8) 



that implies instability if a > — 1. The profile with a < 
— I is stable in the small gap limit. Note that modes with 
q satisfying the resonance condition uja = (or / = 0) 
are marginally stable because T = for them, but V 2 > 
in a neighborhood of the resonance. Therefore, the 
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dependence of T on q should have a two-peak structure 
for any to. As in the case a > 1, the instability occurs 
for any value of e. If a = 1, then we have 



2m 



± 



4TO 2 



(to 2 + g 2 ) 2 



4/x 



(9) 



In this case, the dependence T 2 (q) also has a two-peak 
structure because T ~ at the resonance but T 2 > in 
its neighborhood. The instability is always present for 
any finite value of e. 

Our explicit solution shows that, if a > —1, the insta- 
bility always occur for disturbances with q and to close to 
the condition of magnetic resonance, lua — 0. The axial 
field cannot suppress the instability which occurs even if 
B z is significantly greater than B v q. 



B. Numerical results 
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FIG. 1: The growth rate as a function of q for e = 0.1 and 
a — 1. The panels correspond to m = 1 and m = 100. The 
horizontal axis has different scales in the panels. 



In spite of the various approximations which have been 
done, the picture emerging in the previous session gives 
a qualitatively correct account of the MHD spectrum. In 
order to show this we solved numerically Eq. (5), assum- 



ing a = 1 so that 



The results for other profiles 



of B v are qualitatively similar. Eq. (5) together with the 
boundary conditions is a two-point boundary value prob- 
lem which can be solved by using the "shooting" method 
19]. In order to solve Eq. (5), we used a fifth-order 
Runge-Kutta integrator embedded in a globally conver- 
gent Newton-Rawson iterator. We have checked that the 
eigenvalue was always the fundamental one, as the cor- 
responding eigenfunction had no zero except that at the 
boundaries. 

Fig. 1 exhibits the growth rate of instability as a func- 
tion of q in the case when the toroidal field is stronger 
than the axial one (e = 0.1). We plot T for two values of 
the azimuthal wavenumber, to = 1 and m = 100. Calcu- 
lations confirm that only the modes are unstable with the 
axial wavevectors q close to the condition of the magnetic 
resonance. The resonance values of q = —m/e are 10 and 
1000 for to = 1 and to = 100, respectively. Also, in com- 
plete agreement with the analytic results (see Eq. (9)), 
the growth rate goes to at the resonance but T 2 can be 
positive in its neighborhood. The dependence in Fig, 1 is 
very sharp: the ratio S of the half-thickness of the peak 
to q = —m/e, corresponding to the resonance, is ~ 2 
for to = 1 but rapidly decreases and reaches ~ 0.02 for 
to = 100. The maximum growth rate slowly increases 
with to and becomes ~ 1 for large m that corresponds to 
the growth time of the order of the Alfven crossing time. 
In Fig. 2, we plot the dependence of T 2 on q for the same 
a and e = 10. Qualitatively, the behavior of T 2 is similar 
to that shown in Fig. 1: only modes with q close to the 
magnetic resonance can be unstable, the corresponding 
range of q is narrow and the instability has a resonance 
character, a two-peak structure of T 2 near the resonance, 
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FIG. 2: Same as in Fig. 1, but for e 
spond to m — 10 and m — 200. 
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10. The panels corre- 



the maximum growth rate increases with to, etc. Numer- 
ically, however, the results differ substantially. The reso- 
nance peaks are much sharper for e = 10. For example, 
S is ~ 0.2% and ~ 0.1% for to = 200. The maximum 
growth rate is approximately 10 times lower than in the 
previous figure but still is sufficiently high. Note that, 
generally, disturbances with such small wavelengths in 
the y-and z-directions can be influenced by dissipation 
(viscosity, resistivity). In astrophysical bodies, however, 
the ordinary and magnetic Reynolds numbers are huge 
and even disturbances with to ~ 10 2 — 10 can be treated, 
neglecting dissipation. 
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FIG. 3: Evolution of the mean kinetic density as a function 
of the Alfven travel time in the azimuthal directions for all 
the three models, e = (solid line), 1 (dashed), and 2 (dot- 
dashed) 
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FIG. 4: Radial profile of the eigenfunctions for the fastest 
growing modes excited in the simulations from the linear anal- 
ysis. Solid line represent the e = model with m = 1 and 
q ~ 20. Dashed line and dot-dashed line represent the e = 1 
case for m = 4 and m = 6 respectively at the resonance. 



IV. DIRECT NUMERICAL SIMULATIONS 

It is not difficult to realize this type of instability in 
numerical simulations, at least for moderate values of to. 
In particular, we solved the ideal MHD simulation by 
means of the ZEUSMP code [20[ in the limit of subther- 
mal field. Our setup consists in an isothermal cylinder 
with a radial extent from Sj n to s ou t and vertical size h 
and solve the time dependent ideal MHD equations with 
periodic boundary conditions in z, reflection in s and pe- 
riodic in ip and a resolution ranging from 120 3 to 240 3 
all the directions. The azimuthal field in the basic state 
is taken in the form 

B v = b (s/s ) exp[-(s - s ) 2 /d 2 } (10) 

with bo being a normalization constant; the axial field 
B z is a constant whose value can be varied. In the ba- 
sic state, the Lorentz force is balanced with a gradient 
of pressure, and we have checked that our setup was nu- 
merically stable if no perturbation was introduced in the 
system. For actual calculation we have chosen h = 10, 
Sin = 1.5, s out = 3, so = 2 and d 2 = 0.15; the sound 
speed is assumed to be much larger than the Alfen speed 
(« ten times), in order to compare our results with the 
linear analysis of the previous session obtained for an in- 
compressible plasma. After few time steps we perturbed 
the density with random perturbations in order to excite 
the unstable modes and study their evolution. In the 
case of e = the spectrum is dominated by the m = 1 
mode during the linear phase and we obtain T w 11.7 
for the growth rate in units of the Alfven travel time in 
the azimuthal direction. In order to compare this value 
with the the linear spectrum we explicitly solved Eq.Q 
for our basic state (fT0|) for various values of to and q 
obtaining T m 13.5 for the fastest growing modes for 
the vertical wave numbers excited in the numerical sim- 
ulations according to the spectral analysis. We found 
about ~ 15% difference with the linear result, we think 
this discrepancy is acceptable as 3D simulations are usual 
rather diffusive and one expects that the actual growth 
rate should be smaller than the one obtained from linear 
analysis. Similar considerations apply for the e > 1 cases. 
For instance for e = 1 we find the the fastest growing 
mode has V « 1.54 with to = 4 and m = 6 both excited, 
while the growth rate obtained from the linear analysis 
predicts T s» 1.45. The model with e = 2 has instead 
rn = 9 as the fastest growing modes and also in this case 
the difference with the linear analysis is about 10 — 15%. 
The eigenfunctions corresponding to the fastest growing 
modes for e = and e = 1 are depicted in Fig. (01. In 
Fig- (|3j) the evolution of the mean kinetic energy are plot- 
ted as a function of the Alfven travel time. The solid 
line is for e — 0, while the dashed is for e = 1 and the 
dot-dashed for e = 2. Note that E ax /E tor ~ 13 for model 
e = 1 and E ax /E t0 r ~ 42 for model e = 2 in our setup. 
The growth time for model e = is of the order of the 
Alfven crossing time, while it is significantly longer for 
models e — 1 and e = 2. Nevertheless, the key point that 
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FIG. 5: The density for model e = 2 during the unstable evo- 
lution, around 4a = 7.3, along the plane z = as a function 
of radial and azimuthal coordinate. The presence of a higher 
m mode around m ~ 9 is clearly visible. The domain along <p 
is 2-7T and the resolution of the simulation along the (z, r, tp) 
box is 240 x 120 2 . 



should be stressed here is that the strength of the (turbu- 
lent) magnetic energy and turbulent kinetic energy at the 
beginning of the non-linear phase is essentially the same 
for all the three models. Moreover, in the presence of a 
nonzero axial field the corresponding spectrum along the 
vertical direction shows a specific excited mode, so that 
the resonance condition q ~ —m/e is satisfied. For model 
e — 2 for instance, q w 4 — 5, for the radial component 
of the magnetic field during the linear evolution. Fig.© 
shows the occurrence of high m modes for the density in 
the (s, (p) plane for e = 2 for a 240 x 120 2 simulation. It 
is difficult to reproduce the instability for much higher 
values of e. As it is clear from Fig.© the width of the 
resonance is quite narrow in this case, the growth rate is 
significantly different from zero only for very large values 



of m and the resolution in all three directions needed to 
reproduce the instability can be extremely large. 



V. ASTROPHYSICAL IMPLICATIONS AND 
CONCLUSIONS 

In this paper, we revisited the stability properties of 
the screw pinch, a problem which has received consider- 
able attention in the past in the context of MHD plasma 
stability for thermonuclear fusion. As it was pointed 
out by Freidberg [17| . Eq.(2) describes various types of 
modes which can become unstable under certain condi- 
tions. The basic properties of the unstable modes are 
similar to those of quasi-kinks and quasi-interchanges ob- 
tained by [U [22[ for compressible plasma. However, as- 
trophysical condition like those of stellar interior imply 
a high /3 plasma parameter, a regime which is very far 
from the typical laboratory conditions. To the best of 
our knowledge, an instability of this type has not yet 
been extensively studied for a pressure balanced mixed 
poloidal/toroidal field configuration in the incompress- 
ible limit, an approximation which can be applied to var- 
ious astrophysical situations. The following properties 
characterize the instability in this case: i) the instability 
does not occur for a current-free magnetic configuration; 
ii) it can arise on a time scale comparable to the Alfven 



time scale whereas the growth rate calculated by [22 1 
is an order of magnitude lower, at least; iii) the eigen- 
functions for high values of m have a resonant character 
being very localized as shown in Fig.Q for m = 6; iv) 
the dependence of the growth rate on m seems also to be 
rather peculiar. In the case of the instability described 
in [22] , unfortunately, the growth rate is calculated only 
in the so called tokamak approximation B v /xB z -C 1 
(see Eqs.(30)-(31) by [22j]) and increases approximately 
proportional to m or even faster. In our case, the depen- 
dence on m is qualitatively different because the growth 
rate saturates with m very rapidly, as noticed in the nu- 
merical investigation and in the approximate expression 

©■ 

In spite of these differences, quasi-kink and quasi- 
interchange instabilities obtained by [U [22[ also have 
the typical double-peak structure depicted in Fig.(l) and 
Fig. (2) as a function of the the axial wavevector. 

Note that the basic state in our model is characterized 
by the negative pressure gradient in some fraction of the 
volume, at least. Indeed, hydrostatic equilibrium with 
the toroidal field (10) implies that 



dP 
ds 



2ns \ d 2 



(11) 



Then, dP/ds < if d 2 > s(s - s ). The condition 
dP/ds < is required for the development of instabil- 
ity (see, e.g., [23| ). The sign of the pressure gradient is 
important because it determines the destabilizing effect 
in the so called Suydam's criterion [2~H ]. This criterion 
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represents a necessary but local condition for stability 
and it reads in our notations 

sB 2 z (ldh\ 2 dP „ 

if U*J +8 d7 >0 ' < 12 > 

where ft, = sB z /B v is the magnetic shear. In the case 
of the basic state with toroidal field (10), the necessary 
condition for stability is not satisfied in some fraction 
of the volume (for example, in a neighborhood of So). 
This violation of the stability condition (TT2"j) is actually 
indicating the presence of at least some unstable mode 
in the system. 

Stability properties of magnetic configurations are of 
great importance for various astrophysical applications. 
For instance, it is widely believed that magnetic fields 
play an important role in the formation and propaga- 
tion of astrophysical jets providing an efficient mech- 
anism of collimation through magnetic tension forces 
(e.g., 25]). Polarization observations provide informa- 
tion on the orientation and degree of order of the mag- 
netic field in jets. It appears that many jets can develop 
relatively highly organized magnetic structures. To ex- 
plain the observational data, various simplified models 
of three-dimensional magnetic structures have been pro- 
posed. Typically, the magnetic field can have both lon- 
gitudinal component and substantial toroidal component 
in the core region (see, e.g., [2a|). The mechanisms re- 
sponsible for generation of the magnetic field in jets are 
still unclear. Since the origin of jets is probably relevant 
to MHD-processes in magnetized plasma, their magnetic 
fields could be generated during the process of jet forma- 
tion (see, e.g., [27j ) or, alternatively, it can be generated 
by the dynamo mechanism [28| when the jet propagates 
in the interstellar medium. In both cases, the stability is 



a crucial issue for the properties of the jet. For instance, 
the origin of relatively small scale structures within the 
jet can be attributed to different instabilities arising in 
jets, including the one considered in our study. Magnetic 
structures that appears as a result of the development 
of instabilities can manifest themselves in polarization 
observations of the jets. 

The considered instability can play an important role 
in magnetic stars where it can affect the magnetic field 
in stably stratified regions. Spruit [2!| reviewed various 
types of instabilities that are likely to intervene in a mag- 
netized radiative regions of stars, and he concluded that 
the strongest among them are those which are related to 
the instability of magnetic configurations. According to 
p9j . turbulence generated by such instability can drive 
a genuine dynamo in stellar radiative zones (see, how- 
ever, [30]). Understanding the conditions required for 
the instability is, therefore, crucial for dynamo models in 
stably stratified zones of stars. 

This type of magnetic instabilities can be of interest 
also for neutron stars where the magnetic field reaches 
an extremely high value ~ 10 13 — 10 14 G. Such a strong 
field can be generated by the turbulent dynamo action 
during the very early stage of evolution (see ^31]) when 
the neutron star is convectively unstable. This unstable 
stage lasts less than ~ 1 min. The further evolution of the 
magnetic field is determined mainly by ohmic dissipation 
but can be affected by current-driven instabilities as well 
[32| because dynamo in the convective zone generates a 
magnetic configuration that is not equilibrium. 
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